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Abstract 

Equilibrium and non-equilibrium growth phenomena, e.g., surface growth, generically 
yields self-affine distributions. Analysis of statistical properties of these distributions ap- 
pears essential in understanding statistical mechanics of underlying phenomena. Here, 
we analyze scaling properties of the cumulative distribution of iso-height loops (i.e., con- 
tour lines) of rough self-affine surfaces in terms of loop area and system size. Inspired 
by the Coulomb gas methods, we find the generating function of the area of the loops. 
Interestingly, we find that, after sorting loops with respect to their perimeters, Zipf-likc 
scaling relations hold for ranked loops. Numerical simulations are also provided in order 
to demonstrate the proposed scaling relations. 
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Self-affine distributions are ubiquitous in many phenomena in nature, such as in growing 
surfaces and interfaces [H El 2] , fractured media [6] , and graphs of two-dimensional tur- 
bulent flows [HE]. Self-affine distributions have also been used as a tool to study scaling prop- 
erties of two-dimensional statistical models by mapping these models to a two-dimensional 
Coulomb gas [8 : , 9j. Moreover, crystal growth, the growth of bacterial colonies, and the for- 
mation of clouds in the upper atmosphere [10] are all examples of non-equilibrium phenomena 
which grow self-affine rough surfaces. The above applications on a fundamental level make 
the surface-growth problem as a paradigm for a broad class of problems in the context of 
non-equilibrium statistical mechanics. 

Self-affine surfaces can be described by their height distribution function. From statistical 
point of view, it is necessary to explore topography of this kind of surfaces. In such surfaces, 
heights are invariant under re-scaling, namely h(r) = b~ H h(br), where H is called the rough- 
ness exponent or the Hurst exponent. It implies that in a self-affine surface, the variance of 
the surface height, i.e., yj ([h(x) — (h)] 2 ), scales as L H , where L is the size of the system and 
average is taken over x. If we require translational and rotational invar iance of the surface 
then the structure function of this surface behaves as 



The above equation gives a simple formula to calculate the roughness exponent. To determine 
that a given surface is self-affine or multi-affine we need to measure the pth. order structure 
function defined by C p (r) = (\h(x) — h(x+r)\ p ) . The exponent hierarchy a p is defined through 
the relation C p (r) ~ r a P H _ The exponent a p varies linearly with p for a self-affine surface. 
For a multi-affine surface, instead, it would vary non-linearly with p |llj . The Fourier space 
counterpart of the structure function is Fourier power spectrum S(q) = (|/i(q)| 2 ), where 
/i(q) is the Fourier transform of h(r). Equation (jl.ip gives the scaling relation for the power 
spectrum, i.e., <5(q) ~ Iql -2 ^ 1 "^-*, for small values of q or large values of r. One way to 
generate a Gaussian ensemble of self-affine surfaces is by taking each Fourier height as an 
independent Gaussian random variable with variance given by S(q) ~ \q\~ 2 ( 1+H \ In other 
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words, 

k 



P{h} ~ exp 



d 2 qq 2(1+ ^q^q 



ri.2 



where A = 1/a is the high-momentum cutoff and k is the stiffness. A family of self-affine 
surfaces having all the required properties can be generated by the above distribution. For 
rough surfaces with unbounded heights we have < H < 1, where the higher H is related to 
smoother surface with hills. In a self-affine distribution H > 1/2 (H < 1/2), it implies positive 
(negative) correlations among the increments of the generated values, H = 1/2 means that 
the statistics of the surface follows that of a Brownian motion. At H = 0, it is possible to write 



Eq. (|1.2p in the real space by using ordinary derivative P{h} ~ exp — -| J Q L d 2 x((9 x ^) 



For 



the general case we should replace ordinary derivative with the fractional one, that is, P{h} 



exp 



| J L d 2 x.h(—V 2 ) 1+H h , where the fractional derivative is defined by (— V 2 ) 1+H e* q - X = 
_|^|2+2_ff e iq.x j-£ or more details see [12]). 

The contour lines that are generated by a cut through the surface at a certain height are 
important in characterizing self-affine surfaces. 

In Fig. 1 we plotted an example of the set of contour lines of a rough surface. The 
statistical properties of contour lines of rough surfaces show fractal behavior. The accepted 
fractal dimension of a contour line D = ^-rr- was found by Kondev and Henley (KH) by 
using scaling arguments [13]. Recently, Schramm and Sheffield |14j proved rigorously that 
the contour lines of Gaussian free field with H = are conformally invariant with fractal 
dimension D = |, which is in agreement with the KH result. Conformally invariant curves in 
statistical physics can be investigated by Coulomb gas field theory [8]. The most well-known 
loop model that can be investigated by this field theory is the O(n) model, which can be 
defined on the honeycomb lattice as follows: take the ensemble of loops on the honeycomb 
lattice so that the generating function of the model is given by Z = J^x'ra^, where 
and / are number of the loops and bonds, respectively, n is the weight of each loop, and 
x is the weight of each bond. At the critical point, this loop model can be investigated, 
after mapping the loop model to the solid on solid (SOS) model [8], by the free field theory 



Figure 1: Small part of contour lines of a rough surface with size 3000 2 and H = 0.5; by 
zooming in on the picture one can see many small loops. 

P{h} ~ exp[— | Jq d 2 x(dh) 2 ]. It is also possible to find the scaling exponents of conformal 
curves by the above field theory [8]. Since the height ensemble of a rough surface is not 
conformally invariant, rigorous investigating of their contour lines is more difficult than the 
Coulomb gas case. Indeed, one can not employ the powerful tools of conformal field theory 
(CFT) to study this system. For a rough surface with a generic H there is no rigorous proof 
for results obtained by KH [15]. Nonetheless, it seems that the contour line ensemble shows 
scaling properties similar to the conformal curves encountered in some models such as the 
contour lines of tungsten oxide (WO3) [16] and KPZ surfaces |17j . 

In this paper, by using techniques which are common in the realm of Coulomb gas field 
theory, we introduce new scaling laws for some properties of contour lines of self-affine rough 
surfaces. The scaling properties of the cumulative distribution of the number of contours 
versus the area of the contours and the size of the system are also obtained. In addition, we 
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find a close relation between the cumulants of A, the area of contour lines, and the eigenvalues 
of the fractional Laplacian. Finally, we introduce the scaling property of ranked contour lines 
versus both rank and system size (Zipf's law). Numerical simulations are also provided to 
substantiate our analysis. 

2 Numerical methods 

To generate self-affine rough surfaces in our numerical simulations, 

we have used the successive random addition method [H]. In our simulations we have 
generated surfaces of size LxL with L G {400, 600, 800, 1200, 2000, 3000, 4000}. To investigate 
the effect of roughness exponent on the scaling relations we used several values of H £ 
{0.3, 0.4, 0.5, 0.6, 0.7}. In each case, all calculations have been averaged over 200 realizations. 

To generate the loops in the contour lines we used a contouring algorithm that treats 
the input matrix as a regularly spaced grid. The algorithm scans this matrix and compares 
the values of each block of four neighboring elements (i.e., a plaquette) in the matrix to the 
contour level values. If a contour level falls within a cell, the algorithm performs a linear 
interpolation to locate the point at which the contour crosses the edges of the cell. The 
algorithm connects these points to produce a segment of a contour line. After generating 
the contours of a given surface, in order to eliminate the effect of the edges of the lattice 
we have excluded the contours crossing the edges of the lattice. To show the goodness of 
the fits and consistency of our simulations with theory, we used the following three different 
methods for estimating the exponents: (a) we numerically calculated local slops of the curves 
by fourth-order numerical differentiation for non-uniform data points; e.g., in the case of 
Eq. (|3.2p . derivation of log 10 (A^4) relative to log 10 A (b) We present some of the curves 
by dividing both sides of a scaling relation to the claimed power law to show how seriously 
they are aligned or how they deviate from a horizontal line, e.g., Fig. 2. And, (c) We used 
Bayesian analysis without prior distribution, namely Likelihood analysis [191 1201 l2"Tj [22] to 
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Figure 2: Scaling relation for cumulative distribution of areas. Curves show J^-a/i f° r rou gh 
surfaces with size 4000 2 and different roughness exponents. 



calculate the accuracy of the exponent generated from our numerical results. 



3 Cumulative distribution of area 

A key difference between the contour lines in Coulomb gas field theory and the self-afnne 
rough surfaces is in the fractal dimension of the set of all contour lines. For a given self-afnne 
rough surface, this fractal dimension is d = 2 — H. It is well-known that many of the scaling 
relations in Coulomb gas field theory remain unchanged just by substituting this d as the 
dimension of our set. To give an example, let us define the fractal dimension of a contour 
line D as I ~ R D , where I is the the perimeter of the contour and R is the radius of gyration. 
Moreover, the probability of finding a contour loop with length lis Nr ~ l~ T . One can show 
that there is a hyperscaling relation between the scaling exponents D and r as follows: 

D(r-l)=d, (3.1) 

which is exactly the same as the hyperscaling relation for domain walls in statistical models 
[23]. Following KH, the cumulative distribution of the number of contours with area greater 
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than A has the following form: 

N * ~ ^5- (3-2) 

This gives the right answer for Coulomb gas loops with zero roughness exponent [23]. In the 
rest of the paper using new conjectures we will demonstrate some other evidences to support 
the above relation. This in turn leads us to several new scaling relations. 

We checked Eq. ()3.2p by using numerical simulations for different H's, see Fig 2. As 
is shown, we plot logio( ^d/2 ) versus log 10 (-A) to show how seriously they follow Eq. (|3.2|) . 
The straight horizontal curves exhibit that the proposed scaling relation is preserved up to 
2 orders of magnitude of A. As is seen, in the case of H = 0.3 we have a small deviation 
from the proposed exponent at large values of A, which is related to finite-size effects. For 
a given lattice size and for small values of H, there are not so many large contour lines, but 
we have many small ones. This is led by the nature of self-affinity at small Hurst exponents. 
There are no deviations when we increase H (Fig. 2). In Table [IJ we report the best fit values 
calculated by the likelihood analysis pH [201 EE1 E2] at 68.3% and 95.4% confidence levels. 



H 


Theory 


Local exponent (lex) 


Local exponent (2a) 


0.3 


-0.850 


-0.840 ±0.010 


-0.840 ± 0.020 


0.4 


-0.80 


-0.795 ± 0.006 


-0.795 ± 0.020 


0.5 


-0.750 


-0.752 ±0.010 


-0.752 ±0.025 


0.6 


-0.700 


-0.703 ± 0.005 


-0.703 ± 0.020 


0.7 


-0.650 


-0.652 ± 0.005 


-0.652 ±0.020 



Table 1: The best fit values of exponent d/2 derived by using the likelihood method, a 
denotes standard deviation of each calculated exponent. 

For loops corresponding to surfaces with H = 0, using Coulomb gas techniques, Cardy 
and Ziff showed that C has the universal form as a function of the system size L for different 
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critical statistical physics models [24J. To calculate C, Cardy and Ziff evaluated the total 
area inside all loops using two different methods, and then they found the universal form of 
C. Inspired by this method, we argue to give some new scaling relations for contour lines of 
self-affme rough surfaces. Using Eq. (13, 2p it is straightforward to show that (Atot) = CL H , 
for < H < 1, and for H = it has a logarithmic form. 

Let us consider a typical point x with height h above the horizon (we cut our self-affine 
surface by a plane). If we draw a circle of radius h l l H around the point, since we are dealing 
with a rough surface, all points inside the circle will be above the horizon. In other words, 
inside the loop is a compact region with the fractal dimension 2. Since the fractal dimension 
of the clusters is 2, one could obtain the total area of the clusters proportional to the area 
of the system. This is just a lower bound for the interior areas of the loops — see [25]. In 
addition, it is also possible to see from simulation that by cutting the surface from the average 
height one could get always clusters of the order of the system size. Thus one can get the 
following scaling relation for C with respect to the system size: 



This indicates that the number of contour lines with area greater than A per total length of 
all contours, i.e., L 2 ~ H , is independent of the system size. It is also worth noting that the 
cumulative distribution of the contours with area A is independent of a, the ultraviolet cutoff, 
which is another length scale. Our simulations confirm the validity of the scaling relation 
(|3.3p for different values of H, see Fig. 3. 

The above result is also useful to get another nontrivial equation for contour lines. To 
calculate the total area we can use the formula 



in which d(r) is the minimum number of loops which must be crossed to connect r to the 
edge of the lattice. Since the total area inside the loops is proportional to the area of the 



C ~ L 



2-H 



(3.3) 




(3.4) 
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Figure 3: Scaling relation for the coefficient of cumulative distribution of areas for rough 
surfaces with respect to system size for different values of H as shown in the graph. Slops of 
the curves from top to bottom are given by 1.66 ± 0.04, 1.60 ± 0.03, 1.53 ± 0.03, 1.46 ± 0.04, 
and 1.37 ± 0.04. 

system, we conclude 



This is reminiscent of the height correlation function in the self-affine rough surfaces. For 
H = the relation is logarithmic and was proved explicitly in [24]. These results show 
that one may investigate contour loops of rough surfaces by defining currents for the loops. 
Again, by analogy with the Coulomb gas methods one can define J^(x,y) ~ ^/g-^jTr— as the 
current density of loops. This is a natural candidate if we imagine that the height function 
is extended to the two-dimensional manifold in such a way that it is constant within each 
plaquette. Normalization with respect to width is necessary because we have a rough surface 
where width is changing by size. This definition for the current density means that we can 
map our height model to the contour lines or vice versa. Since iso-height lines have the same 
role as the domain walls between the positive and negative heights, the directional derivative 
of h along a contour must be zero and it must vary along a line normal to the contour. Using 




(3.5) 
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the above function to parameterize the geometry of contour line, it is possible to write 

A = ~\j J\ x ~ X '\KV - i/Vvfo y)Jyi x> i y')dvdv'. (3.6) 

This equation is independent of our definition of currents. By using simple dimensional 
analysis, it is not difficult to find our special normalization, i.e., pj-. One can check that 
Eq. (|3.6p gives < A to t >~ L 2 . Using Eq. (13. 6h . inspired by Cardy's argument [23], we find 
the generating function of the cumulants of area of contour loops. The argument for getting 
cumulants of area is as follows. For the simplicity, we use the Dirichlet boundary condition, 
h =const, on the boundary of the system, which means that loops do not cross the boundary. 
After integration by parts, Eq. (|3.6p gives A ~ J d 2 rh 2 (r). In simulation and experiment 
there are many curves emerging from the boundary and going back to another point in the 
boundary; therefore, there will be no exact Dirichlet boundary condition. 

However, as we will see in the simulations, many of our scaling relations, especially the 
distribution of contours, are independent of the boundary conditions. By using the real space 
representation of the height distribution and the Gaussian integral, one can derive 

A e \(-ff+^H , 2qu wl/2 

{ ' det (_02+3/r)-i/2 ■ V' 1 ' 

where k is the stiffness and u is an auxiliary field. One can write the above equation as an 
infinite sum by using Fourier transform 



2guL< 



(e~ uA ) = e 2Z -»"V fcA - 7, (3.8) 

where are the eigenvalues of the fractional Laplacian with Dirichlet boundary condi- 

tions [12] . Expanding Eq. f|3.8j) gives the higher cumulants of A, 

(An~L 2 ?J2^r- (3-9) 

The sum is convergent for all values of p and H except p — 1 = H = 0, which is logarithmic 
with respect to L. To check the above equation we calculated (A p ) for surfaces with different 
roughness exponents and different sizes. For p = 1 all of the surfaces have (A 1 ) ~ L e with 

10 




Figure 4: Top: Moments (A p ) of loop areas of rough surfaces versus system size; for p = 1 
and 5. Here we have the exponents 6>i = 2.02 ± 0.06, <9 3 = 6.05 ± 0.55, and <9 5 = 11.03 ± 1. 
Bottom: 9 p versus p for surfaces with size 4000 2 with different roughness exponents H 
indicated on the graph. 
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9 = 2 ±0.05. For higher moments one can write (A p ) ~ L p with 9 P ~ 9\p. For surfaces with 
roughness exponent between 0.3 and 0.7, the exponent 9\ varies from 1.94 to 2.15 . One can 
see in Fig. 4b that all of the 9 p J s are linear with respect to p. The deviation from 9 = 2 could 
be related to our restriction in getting larger sizes in simulation. 

4 Zipf 's law for contour lines 

Another interesting scaling relation is the universality of the distribution of the ranked 
loop perimeters, which is named Zipf's law [26]. Following [26} I27j. the average perimeter 
of the nth largest cluster can be found by Eq. (|3.2p . which is called by Mandelbrot the Zipf 
distribution 

T D 

k ~ — , (4.1) 

n d 

where d is the fractal dimension of all loops and D is the fractal dimension of one of loops. 
We should emphasize that we have normalized the equation with the appropriate power of 
total number of contour loops, so we ignore here the scaling of the total number of loops |26j . 
We have numerically checked this scaling relation for self-amne surfaces, both with respect to 
rank n and the system size L. As shown in Fig. 5, in three subfigures (for H G {0.3, 0.5, 0.7}) 
for eight different realizations, we presented the log-log plot of -D/d versus n. Here, l n 
shows a scaling relation according to Eq. (|4,ip . For the case of H = 0.3, the scaling relation 
is preserved for over 2 orders of magnitude of n. Since the number of small loops is few, in 
larger values of H = 0.5,0.7, we could see the agreement just for 1 order of magnitude. To 
calculate the exponent in our numerical results let us consider l n ~ n~^. We calculated the 
exponent; Fig. 5 (bottom right) depicts the variation of £ versus H (the average is over 
different realizations). Since in higher Hs we have lower number of loops, thus £ of higher 
values of H have lower accuracy. We also numerically checked the relation of In versus L. In 
the case of H = 0.3, we obtained a scaling relation with exponent 1.38 ± 0.03, which is near 
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H 



Figure 5: Top left (H = 0.3), right (iif = 0.5), and bottom left (H = 0.3): for system size 
of L = 4000, the curves of ranked loop perimeters divided by n- D l d vs rank numbers are 
shown for eight different realizations. Bottom right: the squares stand for the averaged £ 
(the numerical estimated exponent) over 200 realizations. The solid line shows the theoretical 
relation between £ and H, i.e., £ = ^. 
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theoretical value 1.35. In higher values of H, the estimated exponents are not sufficiently 
accurate because the number of loops in smaller system sizes is low. With the same method 
one can find the average area and the radius of gyration as a function of rank 

L 2 L 
A n ~—, Rn~— (4.2) 

n d n<i 

Both of the above formulas are in good agreement with our numerical results. In these 
kinds of scaling relations the error of the estimated exponents for large system sizes are 
considerably small. We believe Eqs (|4,ip and (|4.2p provide a good method to calculate the 
fractal dimension of a single contour as well as the fractal dimension of all contours. 

5 Discussion and Conclusion 

In summary, by using field theory of rough surfaces and considering current for the model, 
we confirmed the previously known scaling relation for cumulative distribution of area. In 
addition, we found a new scaling relation for this distribution with respect to system size. 
Since the action is not translationally invariant and the small momenta are important, nat- 
urally scaling properties depend to the system size. It seems that large momenta do not 
contribute in the scaling properties. Although system is not invariant under homogeneous 
translation, it is not difficult to see that it is invariant under h — ► h + e^a^x 1 ' , which means it 
is inhomogeneously translational invariant. Using inhomogeneous translation one can define 
the currents J^(x,y) r * j c^i/d^ h corresponding to Wilson loops of the theory and re-derive 
the results of Sec. II. Since we only investigated the scaling properties of the contour lines, 
these two different given currents lead to the same scaling relations. 

Considering these currents for contour lines we think that there may be a close relation 
between the statistics of these lines and the eigenvalues of fractional Laplacian. In this 
paper, we discussed leading scaling behavior with respect to the system size, however, to see 
the effect of the eigenvalues of the fractional Laplacian one needs more careful study of the 
amplitudes as well. Since there is no conformal invariance in the height ensemble, finding the 
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exact values of d(r) and C using the techniques of the Coulomb gas is not tractable. 

We confirmed our proposed scaling relations by simulations through cutting a self-affine 
surface at different heights. We have only interpreted the results for the case of cutting the 
surface at its mean height. But we checked also all of the scaling relations for the cases 
of cutting the surface at heights h = {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9}cr, where a is the 
height variance of the surface. We have not seen any meaningful deviation from what we 
obtained for the mean height. 

We also introduced new Zipf-like scaling relations for the contour lines of self-affine rough 
surfaces, and verified them via simulations. We believe the same scaling relations are applied 
to the clusters of rough surfaces but with different exponents. 
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